Cite this Article
Zhang Qing-Yu, Sun Dong-Ke, Zhu Ming-Fang. A multicomponent multiphase lattice Boltzmann model with large liquid–gas density ratios for simulations of wetting phenomena. Chinese Physics B, 2017, 26(8): 084701  
Permissions
A multicomponent multiphase lattice Boltzmann model with large liquid–gas density ratios for simulations of wetting phenomena
Zhang Qing-Yu1, †, Sun Dong-Ke2, Zhu Ming-Fang1
Jiangsu Key Laboratory for Advanced Metallic Materials, School of Materials Science and Engineering, Southeast University, Nanjing 211189 , China
School of Mechanical Engineering, Southeast University, Nanjing 211189 , China

 

† Corresponding author. E-mail: qingyu.zhang1988@foxmail.com

Abstract

A multicomponent multiphase (MCMP) pseudopotential lattice Boltzmann (LB) model with large liquid–gas density ratios is proposed for simulating the wetting phenomena. In the proposed model, two layers of neighboring nodes are adopted to calculate the fluid–fluid cohesion force with higher isotropy order. In addition, the different-time-step method is employed to calculate the processes of particle propagation and collision for the two fluid components with a large pseudo-particle mass contrast. It is found that the spurious current is remarkably reduced by employing the higher isotropy order calculation of the fluid–fluid cohesion force. The maximum spurious current appearing at the phase interfaces is evidently influenced by the magnitudes of fluid–fluid and fluid–solid interaction strengths, but weakly affected by the time step ratio. The density ratio analyses show that the liquid–gas density ratio is dependent on both the fluid–fluid interaction strength and the time step ratio. For the liquid–gas flow simulations without solid phase, the maximum liquid–gas density ratio achieved by the present model is higher than 1000:1. However, the obtainable maximum liquid–gas density ratio in the solid–liquid–gas system is lower. Wetting phenomena of droplets contacting smooth/rough solid surfaces and the dynamic process of liquid movement in a capillary tube are simulated to validate the proposed model in different solid–liquid–gas coexisting systems. It is shown that the simulated intrinsic contact angles of droplets on smooth surfaces are in good agreement with those predicted by the constructed LB formula that is related to Young’s equation. The apparent contact angles of droplets on rough surfaces compare reasonably well with the predictions of Cassie’s law. For the simulation of liquid movement in a capillary tube, the linear relation between the liquid–gas interface position and simulation time is observed, which is identical to the analytical prediction. The simulation results regarding the wetting phenomena of droplets on smooth/rough surfaces and the dynamic process of liquid movement in the capillary tube demonstrate the quantitative capability of the proposed model.

PACS: 47.11.-j;02.70.-c;47.55.-t;68.08.-p
Keyword:multicomponent multiphase lattice Boltzmann model;large density ratio;contact angle;capillary flow
1. Introduction

Over the last three decades, multiphase lattice Boltzmann (LB) models have developed rapidly and have been employed extensively in the field of computational fluid dynamics.[1–8] As a result of LB models not needing to explicitly track interfaces, these models are ideally suited to study physical phenomena in the multiphase flow systems, e.g., bubbles moving in the liquid phase,[9,10] wetting phenomena at the fluid–solid interfaces,[11–15] droplet–wall collisions,[16] and droplet splashing.[17] Among the current multiphase LB models, the pseudopotential LB scheme has attracted a lot of attention due to its computational efficiency and simplicity in coding.[1–5] The pseudopotential LB model imposes non-local interactions among fluid particles. Accordingly, when the interaction potentials are appropriately chosen, different fluid phases can be separated automatically.[8–21] The original multiphase pseudopotential LB model proposed by Shan and Chen[18] is, however, essentially limited to a small density ratio between the two fluid phases, owing to the increased spurious current at the curved interface between different phases.[1–3] As a result, the model becomes inappropriate for some real-world multiphase fluid flows with large density ratios, such as the water-vapor system.

Many efforts have been devoted to improving the performance of pseudopotential LB models, in particular to increasing the liquid–gas density ratio.[21–42] For single component multiphase (SCMP) LB models, there are two typical techniques to enlarge the density ratio of the two fluid phases: (i) introducing realistic equations of state (EOS),[21–23] and (ii) increasing the isotropy order of the interaction force.[24–28] Yuan and Schaefer[21] proposed an LB model in which the pseudopotential equations were derived from the different EOSs, rather than straightforwardly defined as a function of local density. By incorporating the appropriate EOS into the LB models, such as employing the Carnahan–Starling[22] or van der Waals[23] real gas EOS, the liquid–gas density ratio could increase to higher than 1000:1.

As mentioned above, in the original LB models, the spurious current appearing near a curved phase interface deteriorates the stability of calculations and thus restricts the density ratio of the liquid–gas phases to a small range.[24] To solve the problem, researchers have made efforts to seek the best compromise between a large liquid–gas density ratio and small spurious current.[25–28] Shan[25] considered that the spurious currents existing in the multiphase LB models are caused by insufficient isotropy of the discrete gradient operator. He calculated the interaction force with higher order symmetry, and hence the spurious current effectively decreased. Sbragaglia et al.,[26] and Falcucci et al.[27,28] introduced a mid-range potential interaction into the pseudopotential LB models using two layers of neighboring nodes to calculate the interaction force. Using this approach, the spurious currents at the phase interfaces are diminished by the higher isotropy order. Other techniques, aiming at reducing the spurious currents to increase the liquid–gas density ratios of the pseudopotential LB models, include grid refinement,[29,30] adopting proper force schemes,[31–35] and using the multi-relaxation-time scheme.[36–38]

All the above approaches concern SCMP LB models. In contrast, there are few reports on increasing the liquid–gas density ratios of multicomponent multiphase (MCMP) pseudopotential LB models, yet it is well known that many multiphase flows in nature and practice are composed of multiple chemical components, such as water–air and alloy melt–hydrogen gas flow systems. Liu et al.[39] and Chen et al.[40] elevated the liquid–gas density ratio to a value of higher than 100:1 by imposing an interaction force among the same component particles in the MCMP model. Bao and Schaefer[41] also incorporated the interaction force among the same component particles, which was set to be 103 times higher than that among the different component particles in previous MCMP models. Using such an approach, the density ratio was increased to 1000:1. This approach presents a weakness, because the imposed interaction force among the same component particles might give rise to some problems. For example, it affects the solubility of one component in the other components.[1] Moreover, it also limits the achievable contact angle ranges of droplets on solid surfaces.[40] Another important topic in the field of the MCMP pseudopotential LB model is to implement the different molecular weights (pseudo-particle mass) with respect to the different fluid components. The approaches that have been proposed include (i) applying the interpolation technique or modifying the equilibrium distribution function to derive the distribution functions for different fluid components[43,44] and (ii) adopting the finite-different-based LB model.[45] The two approaches are feasible when the time steps of the different fluid components are identical. However, the studies of adopting different time steps with respect to the different fluid components remain limited in the literature.

In the present work, based on the original MCMP pseudopotential LB model developed by Shan and Chen,[18] a modified MCMP pseudopotential LB model is proposed to simulate the liquid–gas flow with large density contrasts. In the present model, the interaction force among the same component particles is not considered. Instead, we employ a scheme that adopts higher isotropy order to calculate the interaction force, and is combined with a different time step method. Using this approach, the liquid–gas density ratio is enlarged, and the spurious current is restricted for the multicomponent multiphase system. The proposed model is used to simulate the wetting behavior of droplets on smooth/rough solid surfaces and the dynamic process of liquid movement in a capillary tube. The simulated intrinsic and apparent contact angles of the droplets on solid surfaces, as well as the simulated time evolution of the moving liquid–gas interface position in a capillary tube, are compared with theoretical predictions.

2. Model description

The techniques used in the present MCMP model involve (i) using two layers of neighboring nodes to calculate the interaction forces between the two components, and (ii) adopting different time steps for the two components.

Figure 1 shows the node distributions for calculating the fluid–fluid cohesion force in the original and present pseudopotential LB models, respectively. As shown in Fig. 1(a), only the nearest nodes (nodes 1, 2, 3, 4) and the next-nearest nodes (nodes 5, 6, 7, 8) are considered in the original model to calculate the cohesion force, yielding the fourth isotropy order (E4).[25] In the present model, however, two layers of neighboring nodes (24 nodes) are used to calculate the cohesion force as shown in Fig. 1(b), obtaining the eighth isotropy order (E8) that can reduce the spurious currents at the phase interface.[25–28,42] In addition, the time step of component 2 (the essential component in the gas phase) is set to be larger than that of component 1 (the essential component in the liquid phase), so that the components with a large pseudo-particle mass contrast can be implemented. This scheme takes into account the fact that the mean free path of gas particles is longer than that of liquid particles. This causes the period of component 2 for accomplishing the processes of propagation and collision to be times longer than that of component 1, where and are the time steps of component 1 and component 2, respectively. Accordingly, the time step ratio ( between the two components is set to be based on the E8 MCMP LB model. It is noted that the computation efficiency of the present MCMP model is inversely proportional to the value of the time step ratio. The pseudo-particle mass ratio between the two components is assumed to be . This magnitude of is selected according to the large weight ratios of fluid particles in real-world liquid–gas systems, such as alloy melt–hydrogen gas flows.

Fig. 1. (color online) Sketch of the node distributions for calculating the fluid–fluid cohesion force in (a) E4 and (b) E8 LB models.

The two fluid components are represented by two sets of distribution functions. The processes of particle propagation and collision are governed by the LB equations using the two-dimension, nine-velocity (D2Q9) scheme:

where , , and are the density distribution function, the time step, and the relaxation time of the σ-th component, respectively. The relaxation times of the two components are adopted to be in the present study. Discrete velocities, , in the D2Q9 scheme are given by
where is the lattice speed, and is the lattice spacing. The equilibrium distribution function, , is expressed as
where the weight coefficients ωi are given as , , and , and is the number density of the σ-th component, which is obtained from . Based on the original MCMP pseudopotential LB model, the time step-related macroscopic velocity is given by
where mσ is the pseudo-particle mass of the σ-th component. In the present study, the time steps of components 1 and 2 are assumed to be and , respectively. Since component 1 accomplishes the propagation and collision processes times with , while component 2 accomplishes it once with , it is considered that the microscopic velocity of component 1 should interact with that of component 2 times. Therefore, in the numerator of the first term on the right-hand side of Eq. (4), the time steps are taken into consideration when calculating the microscopic momentum. in Eq. (4) is the interaction force acting on the σ-th component, which includes the fluid–fluid cohesion force and the fluid–solid adhesion force .

In the original MCMP pseudopotential LB model with E4 isotropy order, the fluid–fluid cohesion force acting on the σ-th component, which generates different fluid phases, is calculated from

where σ and denote the two different fluid components, is the fluid–fluid interaction strength that controls the cohesion force and surface tension between two fluid phases, and and are the pseudopotential functions.

As described above, for reducing the spurious current at the phase interface, the present model adopts the E8 isotropy order to calculate the cohesion force. Accordingly, equation (5) is rewritten as

where and denote the first layer (nodes 1–8 in Fig. 1(b)) and the second layer (nodes 9–24 in Fig. 1(b)) of the neighboring nodes, respectively. The numerical values of the weight coefficients, namely, p1i and p2i, are listed in Table 1.[25–28,42] In the present study, the periodic boundary condition is imposed on the walls of the computational domain. For the domain consisting of computational nodes, the node index is assumed to be in the x-axis direction and , in the y-axis direction, respectively. When the cohesive forces are calculated by using Eq. (6), the neighboring nodes of the fluid particles at nodes (close to the right boundary) involve those located at the positions of , , (the right boundary), and (0, j) (the left boundary), respectively. The same approach is used for the neighboring nodes of the fluid particles located at the layers close to the other boundaries. The pseudopotential functions of the two fluid components are taken as fluid number densities, which are and , respectively.

Table 1.

Weight coefficients of the E8 model.[25–28,42]

.

When the system includes a solid phase, the fluid–solid adhesion force, , acting on the σ-th component is computed from

where is an indicator function that is equal to 1 or 0 for a solid or a fluid neighboring node, respectively. The interaction strength between the fluid and solid phases can be adjusted by the fluid–solid interaction parameter, . In the present work, a bounce-back rule is adopted on the fluid–solid boundaries.

3. Results and discussion
3.1. Spurious current

In the present study, the spurious current is calculated from

The E8 model, having a higher isotropy order than the original E4 model, is adopted to reduce the spurious current that may deteriorate the stability of the LB model. In order to examine the effect of spurious current reduction, the calculations of a stationary droplet equilibrated with gas phase are conducted by using the E4 and E8 MCMP LB models, respectively. For simplicity, the time steps and the pseudo-particle mass values of the two components are set to be = 1 and = 1 in Subsecttion 3.1, respectively. At the start of the simulation, a droplet with an initial radius of 10 lattice units is placed in the center of the computation domain that is composed of 60 × 60 lattice units. After the equilibrium states are achieved, the maximum spurious currents appearing at the liquid–gas interfaces are measured. Figure 2 shows the comparison between the maximum spurious currents obtained from the E4 LB model and those from E8 LB models. It is shown that the maximum spurious currents for both E4 and E8 models increase with the fluid–fluid interaction strength, , increasing. For a given , however, the maximum spurious current in the E8 model is remarkably lower than that in the E4 model. This result accords well with the result given by Shan[25] regarding the SCMP LB model. It is shown that the implementation of the E8 model with higher isotropy order is capable of reducing the spurious current, thus effectively increasing the computation stability of the MCMP LB model.

Fig. 2. (color online) Comparison between the maximum spurious currents varying with the fluid–fluid interaction strength from the E4 model and those from the E8 model.
3.2. Liquid–gas density ratio

In this study, the liquid–gas density ratio is obtained from the calculations of a static droplet surrounded by the gas phase. The computation domain size and the initial droplet radius are identical to those in Fig. 2. After the droplet equilibrates with the surrounding gas phase, the droplet is mainly composed of component 1 with slightly dissolved component 2. In contrast, component 2 is the essential composition of the gas phase with a slight amount of component 1. Therefore, the liquid–gas density ratio is calculated from , where , , , and are the number densities of component 1 and component 2 in the liquid and gas phases, respectively. The liquid–gas density ratio can be affected by several factors, such as the values of the fluid–fluid interaction strength (, the time step ratio (, the relaxation times (τ1 and τ, and the pseudo-particle mass ratio ( of the two components. In the present study, the values of τ1, τ2, m1, and m2 are assumed to be unchanged, while the effects of and are systematically studied.

Figure 3 presents the obtained liquid–gas density ratio () as a function of time step ratio ( with various values of fluid–fluid interaction strength ( by using the proposed model. It is shown that the liquid–gas density ratio increases with increasing when is fixed. Note that the value of controls the separation degree of the two components. As a result, a larger tends to reduce the solubility of component 1 in the gas phase, which accordingly reduces the gas density and elevates the liquid–gas density ratio. On the other hand, for a given , the liquid–gas density ratio increases as becomes larger before it approaches to a relatively stable value. It is found that as the time step ratio increases, the quantity of the fluid particles of component 1 in the gas phase gradually becomes lower, which facilitates the higher liquid–gas density ratio. Under the condition of , the maximum liquid–gas density ratio is achieved. When the time step ratio is assigned as , the quantities of the fluid particles of component 1 in the liquid and gas phases almost remain stable, leading to a nearly unchanged liquid–gas density ratio. In Fig. 3, it is found that the maximum liquid–gas density ratio achieved by the present model is under the conditions of and .

Fig. 3. (color online) Liquid–gas density ratios varying with time step ratio for various fluid–fluid interaction strengths.

As discussed above, the spurious current appearing at the interface is one of the key factors that hinder the liquid–gas density ratio from being elevated. The influences of the fluid–fluid interaction strength and time step ratio on the maximum spurious current are tested, and the results are presented in Fig. 4. The upper-left inset in Fig. 4 illustrates the map of the velocity distribution, indicating that the maximum spurious current appears at the liquid–gas interface. It can be seen from Fig. 4 that the magnitude of the maximum spurious current is strongly dependent on the fluid–fluid interaction strength , but weakly influenced by the time step ratio . Apparently, the maximum spurious current increases with for each . It is found that to obtain stable calculations of a static droplet surrounded by gas phase, the maximum applicable fluid–fluid interaction strength is around 4.6.

Fig. 4. (color online) Plots of maximum spurious current varying with time step ratio for various fluid–fluid interaction strengths.

It should be noted that the maximum liquid–gas density ratios discussed above are for the case of the static multiphase flow. When the external fluid flow is introduced in the computation domain, e.g., the shear flow,[46] the achievable maximum liquid–gas density ratio is found to become lower than that obtained from the static liquid–gas flow.

3.3. Intrinsic contact angles of droplets on a smooth solid surface

The proposed model is used to simulate droplets on a smooth solid surface with various intrinsic contact angles. The simulation domain is composed of 100 × 100 lattice units and the droplet radius is initially taken as 15 lattice units. The time step ratio is fixed as . For a given fluid–fluid interaction strength , the intrinsic contact angle of a droplet on a solid surface is determined by the fluid–solid interaction parameters, namely and . Regarding the contact angle simulations by using the original MCMP LB model reported in the literature, the relationship between and is not unique.[19,44] In the present study, is adopted in the simulations of contact angles, which is identical with the configuration in the work of Huang et al.[47] Figure 5 illustrates droplets on a smooth solid surface with different contact angles simulated by using the corresponding values with a fixed fluid–fluid interaction strength of . The contact angles in Fig. 5 are measured by using the approach proposed by Huang et al.[47]

Fig. 5. (color online) Simulated contact angles with the fluid–fluid interaction strength of and the fluid–solid interaction strengths of (a) , (b) , (c) , (d) , (e) , (f) (red: solid, blue: liquid, green: gas).

Sukop and Thorne[19] and Huang et al.[47] introduced two forms of empirical formulas that are composed of the LB interaction strengths and related to Young’s equation. The empirical formulas can be used to predict the intrinsic contact angles simulated using the original MCMP pseudopotential LB models. Since the present model implements the higher isotropy order and the different-time-step method, the previous empirical LB Young equations in the literature[19,47] are not appropriate for this work. Thus, a new empirical formula of Young’s equation is established based on the simulations of Fig. 5. The interaction strengths, , , and which are related to the interfacial tensions between the fluid–fluid and fluid–solid phases, are rearranged to construct an empirical formula of Young’s equation as follows:

where and are the fluid–solid interaction strengths when the contact angle is 90°. and are the liquid and gas densities, respectively. It is noted that equation (8) has a form similar to the empirical formula suggested by Sukop and Thorne[19] and Huang et al.[47]

Figure 6 shows the comparison between the LB simulations and calculations by Eq. (8) for the contact angles as a function of fluid–solid interaction strength with various values of fluid–fluid interaction strengths . As mentioned above, the magnitude of determines the surface tension and density ratio of liquid–gas phases. It can be seen that the simulated contact angles are in good agreement with those calculated by Eq. (8) for a wide range of contact angles, even though the simulations deviate slightly from the predictions of Eq. (8) in the high contact angle region for the case of , and the maximum difference is around 8.6°. Moreover, as shown in Fig. 6, in the range of , corresponding to the contact angle region of , the influence of on the contact angles is relatively limited. However, as the contact angles increase or decrease, departing from 90° by regulating , the effect of becomes increasingly evident. For example, with fixing , the contact angles remain close to 90°, regardless of the variation of . In contrast, the contact angles can be obviously influenced by with selecting or which produces high or low values of contact angles, respectively. Figure 6 also reveals that the wettability of the smooth surface increases when adopting a larger . This is because a higher magnitude of leads to a stronger adhesion force between the solid surface and component 1 that is the essential composition of the liquid droplet. The stronger adhesion force drives the droplet to spread on the surface and hence a relatively low intrinsic contact angle is obtained.

Fig. 6. (color online) Comparison between LB simulations (symbols) and predictions by the empirical formula of Young’s equation (lines) regarding the droplet contact angles on a smooth surface.

Figure 7 shows the liquid–gas density ratio in the presence of solid phase as a function of the droplet intrinsic contact angle with various values of fluid–fluid interaction strengths . It is seen that for a fixed value of , the liquid–gas density ratio almost remains stable with a slight decreasing trend as the intrinsic contact angle becomes larger. For a fixed intrinsic contact angle, the increase of elevates the liquid–gas density ratio drastically.

Fig. 7. (color online) Liquid–gas density ratios as a function of intrinsic contact angle in the presence of solid phase for various fluid–fluid interaction strengths.

Comparing Fig. 7 with Fig. 3, it is found that the liquid–gas density ratio in the presence of the solid phase is nearly identical to that of the liquid–gas system without the solid phase when the same level of is adopted. For example, when is taken to be 3.0, the liquid–gas density ratios are both around 300:1 for the two systems with and without the presence of a solid phase. However, the applicable value for the solid–liquid–gas system is smaller than that for the liquid–gas system, which hinders the liquid–gas density ratios from being enlarged to a high level in the presence of a solid phase. To explain this phenomenon, the maximum spurious current is analyzed for the solid–liquid–gas system.

Figure 8 illustrates the maximum spurious currents in the presence of a solid phase as a function of the intrinsic contact angle for three different values of fluid–fluid interaction strengths . In the present two-dimensional (2D) study, the maximum spurious current is found to appear at the solid–liquid–gas three-phase contact points as shown in the velocity contour map in Fig. 8. This implies that for the same value of , the maximum spurious current appearing in the solid–liquid–gas system is larger than that in the liquid–gas system. Moreover, it can be seen from Fig. 8 that for a certain , the intrinsic contact angles within produce relatively small values of the spurious current. Note that the intrinsic contact angles in this range correspond to the simulations adopting fluid–solid interaction strengths of . As the contact angles are regulated to higher or lower degrees by increasing and , the spurious current is accordingly elevated, indicating that the spurious current is enlarged by the fluid–solid interaction force. For a certain value of the intrinsic contact angle, the spurious current is apparently enlarged by the increase of . As discussed previously, the local intense spurious currents may lead to non-convergent computations. The magnitude of is thus restricted to a lower level for the solid–liquid–gas system than that for the liquid–gas system, owing to the fact that there is higher maximum spurious current appearing in the solid–liquid–gas system. Since is a dominant parameter responsible for enlarging the liquid–gas density ratios as shown in Figs. 3 and 7, it is evident that the achievable maximum liquid–gas density ratio in the simulations of the solid–liquid–gas system are smaller than that of the liquid–gas system. In the present work, the maximum liquid–gas density ratio obtained from the solid–liquid–gas system is around 600:1.

Fig. 8. (color online) Maximum spurious currents as a function of intrinsic contact angle in the presence of solid phase for three diferent fluid–fluid interaction strengths.
3.4. Apparent contact angles of a droplet on rough solid surfaces

Li et al.[48] performed simulations of the apparent contact angles of droplets on rough surfaces by using two types of SCMP pseudopotential LB models. The intrinsic contact angles of the droplets in their LB models were regulated close to 138°. After the steady states of the droplets were achieved, the apparent contact angles were measured and compared with the predictions of Cassie’s law. To further demonstrate the validity of the present MCMP LB model, the wetting behaviors of droplets on rough surfaces are also simulated. The computation domain consists of 150 × 100 lattice units and an initial droplet with a radius of 30 lattice units is placed on the rough surface. As shown in Fig. 9, the rough surface is composed of square pillars. The width, pillar–pillar spacing, and height of the square pillars are selected as w = 6, s = 3, and h = 20 lattice units, respectively. The initial droplet radius and the geometrical configuration of the rough surface are taken to be identical with those used in the study of Li et al.[48] The time step ratio is set to be that is identical with those used in the simulations of droplets on a smooth surface in Subsection 3.3. With respect to the three liquid–gas regimes determined by , , and , the fluid–solid interaction strengths are regulated as , , and , respectively, to obtain an intrinsic contact angle of 137.7°. After the droplets reach the steady state on the rough surfaces, the apparent contact angles, measured from the simulated droplets shown in Fig. 9, are 142.4°, 140.0°, and 140.0°, respectively. According to Cassie’s law,[49] the apparent contact angle of a droplet sitting on the top of the pillars can be calculated from

where and θ denote the apparent contact angle and the intrinsic contact angle, respectively and f is the area fraction of the liquid–solid contact surface with respect to the horizontally projected plane of the droplet base. In the present (2D) study, according to the shape of droplets on the rough surfaces shown in Fig. 9, f is calculated from . Substituting the area fraction f = 5/7 and the intrinsic contact angle into Eq. (9), the apparent contact angle is computed to be . Thus, the relative errors between the simulation data measured from Fig. 9 and the predictions of Eq. (9) are all less than 3.2%. In the work of Li et al.,[48] the simulated apparent contact angles on the rough surface are all around 146.0°. It can be found that the simulated apparent contact angles in Fig. 9 are also close to the simulation results obtained by Li et al.[48] These obtained agreeable results indicate that the proposed MCMP LB model can reasonably describe the wetting phenomena of droplets on rough surfaces in solid–liquid–gas systems with various interaction strengths.

Fig. 9. (color online) Simulated patterns of liquid droplets sitting on a rough surface. Three cases of simulations have identical intrinsic contact angles of 137.7° by adopting (a) and , (b) and , (c) and (red: solid, blue: liquid, green: gas).
3.5. Dynamics of capillary flow

The above contact angle simulations are all for the static cases. In order to demonstrate the validity of the present MCMP LB model for the simulation of fluid dynamics in the liquid–gas–solid coexisting system, the dynamic process of liquid movement in a capillary tube is simulated. Liquid moves in a capillary tube due to the pressure difference over the curved liquid–gas interface, which is simulated using several types of LB multiphase models, and can be taken as a benchmark for assessing the capability of a multiphase model for the simulation of dynamic problems.[50–52] The position of the moving liquid–gas interface in the capillary tube can be expressed as follows:[52]

where γ is the surface tension, θ is the intrinsic contact angle, r is the width of the capillary tube in the 2D study, ξ is the real-time position of the liquid–gas interface in the capillary tube, and are the dynamic viscosities of the liquid and gas phases, respectively.

In the present study, the computation domain consists of 300 × 25 lattice units. In the middle portion of the computation domain, namely, , placed is a horizontal capillary tube with a width of 17 lattice units. The liquid region is initialized at the position of , and the remaining region is filled with the gas phase. Thus, the initial liquid–gas interface position on the left-hand side, ξ, is x = 100 in the computation domain. The simulation parameters, including the time step ratio and interaction strengths, are assigned to be identical with those of Fig. 9. Consequently, according to the simulations of Fig. 9, the intrinsic contact angles are 137.7° for the three liquid–gas regimes distinguished by , and that determine different surface tensions. Figure 10 displays the time evolutions of liquid movement in the capillary tube when the interaction strengths are taken as and . It is observed that after the calculation is started, the initial plane liquid–gas interface becomes curved, and the liquid moves in the capillary tube gradually towards the right-hand side.

Fig. 10. (color online) Simulated liquid movements in a capillary tube (red: solid, blue: liquid, green: gas).

Since the relaxation times of the two components are adopted as , the dynamic viscosities of the liquid and gas phases are identical, namely, . Therefore, equation (10) can be rewritten for the LB simulations as

where M is a parameter and defined as . As described in Subsection 3.4, with respect to the three liquid–gas regimes determined by , , and , the contact angles are regulated to be identical. Therefore, M in Eq. (11) is a constant in the simulations of the three liquid–gas regimes. It can be seen from Eq. (11) that the variation of the liquid–gas interface position, ξ, with respect to simulation time, t, is linear. The slope of the ξ–t line is proportional to the surface tension, γ. In the present MCMP LB model, the surface tension can be calculated through Laplace’s law by measuring the pressure difference across the liquid–gas interface and the radius reciprocal of the spherical droplets. It is obtained that the surface tensions of the three liquid–gas regimes are , , and for , , and , respectively, showing that the surface tension, γ, increases with the fluid–fluid interaction strength, .

Figure 11 shows the plots of the simulated position of the curved liquid–gas interface in the capillary tube versus simulation time for different values of fluid–fluid interaction strengths, . The symbols in Fig. 11 represent the simulation results, and the lines are the linear fitting plots of the symbols. Note that the linear relation between the liquid–gas interface position and simulation time is observed, which is identical with the analytical prediction of Eq. (11). In addition, it is found that the slope of the ξ–t line increases with the increase of fluid–fluid interaction strength, . The slopes can be determined from the linear fitting plots in Fig. 11, and obtained as , , and for (), (), and (), respectively. Accordingly, the slope of the ξ–t line is approximately proportional to the surface tension, which is identical with the analytical prediction by Eq. (10). Liu et al.[52] performed the simulation of the capillary flow by using a color-fluid multiphase LB model, and also obtained the linear relationship between the liquid–gas interface position, ξ, and time, t. These agreeable results demonstrate the capability of the present MCMP LB model when dealing with the dynamic problems in the liquid–gas–solid systems.

Fig. 11. (color online) Time evolutions of the liquid–gas interface position in a capillary tube for various fluid–fluid interaction strengths (symbols: LB simulations, lines: linear fitting).
4. Conclusions

In the present study, a modified multicomponent multiphase (MCMP) pseudopotential lattice Boltzmann (LB) model with large liquid–gas density ratios is proposed. Two layers of neighboring nodes are employed to calculate the fluid–fluid cohesion force, which can elevate the isotropy order of the interaction force and effectively reduce the spurious current. Different time steps for the two fluid components, with a large pseudo-particle mass contrast, are implemented in the proposed MCMP LB model. It is found that the liquid–gas density ratio can be elevated by increasing at a fixed . For a given , the liquid–gas density ratio increases as becomes larger before it approaches to a relatively stable value. For the liquid–gas flow system without solid phase, the largest liquid–gas density ratio is observed to be under the conditions of and .

The maximum spurious current appears at the liquid–gas interface for the liquid–gas system and at the solid–liquid–gas three-phase contact points for the flow system in the presence of a solid phase. The maximum spurious current is found to be obviously elevated by increasing , but weakly affected by . Moreover, the magnitudes of the fluid–solid interaction strengths, namely and , have less effect on the liquid–gas density ratio, but influence the spurious current clearly. Accordingly, the applicable value for the solid–liquid–gas system is smaller than that for the liquid–gas system, which hinders the liquid–gas density ratio from being enlarged to a high level in the presence of solid phase. As a result, the achievable maximum liquid–gas density ratio is smaller than that obtained from fluid flows in the absence of a solid phase.

Different contact angles of droplets on both smooth and rough solid surfaces are simulated by using the proposed MCMP LB model. The different intrinsic contact angles are obtained by adjusting the interaction strengths. An empirical formula of Young’s equation containing LB interaction strengths is proposed to predict the intrinsic contact angles. The intrinsic contact angles obtained from the simulations and calculated from the suggested LB Young’s equation are in good agreement with each other. With the droplets sitting on the tops of the pillars of a rough surface, the simulated apparent contact angles with various interaction strengths accord well with the predictions of Cassie’s law. Moreover, the dynamic process of liquid movement in a capillary tube is simulated. Under the condition of the liquid and gas phases having an identical dynamic viscosity, the linear relation between the liquid–gas interface position in the capillary tube and simulation time is observed, and the slope of the linear line is approximately proportional to the surface tension. These simulation results of the liquid movement in a capillary tube are consistent with the analytical predictions. The simulations of the static and dynamic problems show that the proposed MCMP LB model is capable of reasonably simulating the wetting phenomena in the solid–liquid–gas coexisting systems with large liquid–gas density ratios.

Reference
[1] Chen L Kang Q Mu Y He Y L Tao W Q 2014 Int. J. Heat Mass Transfer 76 210
[2] Li Q Luo K H Kang Q J He Y L Chen Q Liu Q 2016 Prog. Energy Combust. Sci. 52 62
[3] Succi S 2015 Europhys. Lett. 109 50001
[4] Yong W A Zhao W Luo L S 2016 Phys. Rev. 93 033310
[5] Sun D Zhu M Wang J Sun B 2016 Int. J. Heat Mass Transfer 94 474
[6] Sun D Wang Y Dong A Sun B 2016 Int. J. Heat Mass Transfer 94 306
[7] Sun D Zhu M Pan S Raabe D 2009 Acta Mater. 57 1755
[8] Yang K Guo Z 2016 Phys. Rev. 93 043303
[9] Zheng H W Shu C Chew Y T 2006 J. Comput. Phys. 218 353
[10] Liang H Li Q X Shi B C Chai Z H 2016 Phys. Rev. 93 033113
[11] Zhang Q Y Sun D K Zhang Y F Zhu M F 2014 Langmuir 30 12559
[12] Farokhirad S Morris J F Lee T 2015 Phys. Fluids 27 102102
[13] Patrick Jansen H Sotthewes K Zandvliet H J W Kooij E S 2016 Appl. Surf. Sci. 361 122
[14] Xie C Y Zhang J Y Bertola V Wang M R 2016 J. Colloid Interface Sci. 463 317
[15] Zhang Q Y Sun D K Zhang Y F Zhu M F 2016 Chin. Phys. 25 066401
[16] Wang Y Shu C Yang L M 2015 J. Comput. Phys. 302 41
[17] Raman K A Jaiman R K Lee T S Low H T 2015 Comput. Fluids 107 285
[18] Shan X W Chen H D 1993 Phys. Rev. 47 1815
[19] Sukop M C Thorne D T 2006 lattice Boltzmann Modeling: An Introduction for Geoscientists and Engineers Heidelberg, Berlin, New York Springer 112 115
[20] Khajepor S Wen J Chen B 2015 Phys. Rev. 91 023301
[21] Yuan P Schaefer L 2006 Phys. Fluids 18 042101
[22] Chen X P Zhong C W Yuan X L 2011 Comput. Math. Appl. 61 3577
[23] Hu A Li L Chen S Liao Q Zeng J 2013 Int. J. Heat Mass Transfer 67 159
[24] Connington K Lee T 2012 J. Mech. Sci. Technol. 26 3857
[25] Shan X 2006 Phys. Rev. 73 047701
[26] Sbragaglia M Benzi R Biferale L Succi S Sugiyama K Toschi F 2007 Phys. Rev. 75 026702
[27] Falcucci G Chibbaro S Succi S Shan X Chen H 2008 Europhys. Lett. 82 24005
[28] Falcucci G Ubertini S Succi S 2010 Soft Matter 6 4357
[29] Yu Z Fan L S 2009 J. Comput. Phys. 228 6456
[30] Yu Z Yang H Fan L S 2011 Chem. Eng. Sci. 66 3441
[31] Huang H B Krafczyk M Lu X Y 2011 Phys. Rev. 84 046710
[32] Sun K Wang T Y Jia M Xiao G 2012 Physica 391 3895
[33] Lycett-Brown D Luo K H 2015 Phys. Rev. 91 023305
[34] Zarghami A Looije N Van den Akker H 2015 Phys. Rev. 92 023307
[35] Hu A Li L Uddin R 2015 J. Comput. Phys. 294 78
[36] Li Q Luo K H Li X J 2013 Phys. Rev. 87 053301
[37] Zhang D Papadikis K Gu S 2014 Int. J. Multiphase Flow 64 11
[38] Xu A Zhao T S An L Shi L 2015 Int. J. Heat Fluid Flow 56 261
[39] Liu M Yu Z Wang T Wang J Fan L S 2010 Chem. Eng. Sci. 65 5615
[40] Chen L Kang Q Tang Q Robinson B A He Y L Tao W Q 2015 Int. J. Heat Mass Transfer 85 935
[41] Bao J Schaefer L 2013 Appl. Math. Model. 37 1860
[42] Yang J H Boek E S 2013 Comput. Math. Appl. 65 882
[43] Mccracken M E Abraham J 2005 Phys. Rev. 71 046704
[44] Chai Z H Zhao T S 2012 Acta Mech. Sin. 28 983
[45] Zheng L Guo Z Shi B Zheng C 2010 Phys. Rev. 81 016706
[46] Wang N Liu H Zhang C 2016 J. Comput. Sci. 17 340
[47] Huang H B Thorne D T Schaap M G Sukop M C 2007 Phys. Rev. 76 066701
[48] Li Q Luo K H Kang Q J Chen Q 2014 Phys. Rev. 90 053301
[49] Frank X Perré P Li H Z 2015 Phys. Rev. 91 052405
[50] Pooley C M Kusumaatmaja H Yeomans J M 2009 Eur. Phys. J. Spec. Top. 171 63
[51] Diotallevi F Biferale L Chibbaro S Lamura A Pontrelli G Sbragaglia M Succi S Toschi F 2009 Eur. Phys. J. Spec. Top. 166 111
[52] Liu H Ju Y Wang N Xi G Zhang Y 2015 Phys. Rev. 92 033306